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CHAPTER  I 
INTRODUCTION 

There  has  been  considerable  concern  in  the  past  few  years 
about  the  formulation  of  efficient  and  accurate  numerical  schemes 
suitable  for  integrating  numerical  weather  prediction  models  on  the 
globe.  This  interest  has  grown  not  only  from  the  need  to  have  weather 
predictions  for  the  entire  globe  but  also  from  the  realisation  that  for 
forecast's  of  more  than  a few' days  the  introduction  of  artificial  bound- 
aries anywhere  will  deteriorate  the  forecast  product. 

The  introduction  of  standard  numerical  techniques  into  the 
global  prediction  problem  has  proven  to  be  far  from  straightforward. 
The  difficulty  is  related  to  the  singularity  of  the  spherical  coordinate 
system  at  the  poles.  The  precise  singularities  at  the  poles  can  be 
avoided  e.g.  Merilees  (1973)  but  a basic  difficulty  remains  , namely  , 
the  convergence  of  the  meridians  requires  excessively  small  time 
steps  to  ensure  the  numerical  stability  of  the  model.  This  requirement 
is  costly  in  terms  of  computer  time. 

A second  , more  general  , concern  has  been  with  the  accu- 
racy of  numerical  estimation  of  derivatives  . Numerous  studies  have 
shown  that  the  standard  second  order  approximations  to  derivatives 
on  standard  gridlengths  are  simply  not  sufficiently  accurate  to  make 
numerical  error  a secondary  contributor  to  forecast  error.  For  exa- 
mple Chouinard  and  Robert(1972)  have  shown  that  with  400km  grid  and 
a second  order  finite  difference  scheme  , these  errors  may  account 
for  20%  of  the  rms  geopotential  error  in  a 3 6- hour  forecast. 


For  both  these  reasons  there  has  been  considerable  interest 
in  the  development  of  global  spectral  models  in  terms  of  spherical 
harmonics.  Such  a numerical  basis  has  the  double  advantage  of  a lack 
of  singularity  in  the  numerical  basis  and  no  linear  truncation  error 
in  the  computation  of  derivatives.  Further  , by  use  of  numerical  rather 
than  analytical  transforms  it  is  possible  to  implement  such  models 
with  reasonable  efficiency  , certainly  on  3rd  generation  computer 
systems  presently  available.  However  such  a numerical  basis  does 
require  a considerably  amount  of  computation  per  degree  of  freedom 
and  therefore  it  is  wise  to  continue  to  investigate  the  viability  of 
alternative  schemes  which  are  less  costly  than  the  complete  spectral 
method  yet  have  at  least  some  of  their  advantages. 

These  consideration  led  Orszag(l972)  to  discuss  the  use  of 
the  linear  property  of  the  spectral  method  without  going  to  a complete 
spectral  model  to  get  a method  with  considerable  operational  advant- 
ages over  the  spectral  approximation.  This  method  , known  as  the 
pseudospectral  method  , retains  the  traditional  grid  points  for  the 
representation  of  the  meteorological  fields  , but  expresses  these  fields 
as  finite  Fourier  series  for  the  purpose  of  estimation  of  derivatives. 

As  such  , the  method  can  take  advantage  of  the  fast  Fourier  transform 
(FFT)  to  estimate  derivatives.  Orszag  in  the  above  mentioned  paper 
shows  that  the  error  of  the  pseudospectral  approximation  compared 
to  the  spectral  (Galerkin)  approximation  in  some  simple  models  are 
similar  , despite  the  inclusion  of  aliasing  terms  in  the  pseudospectral 
approximation. 


- 3 - 

As  mentioned  , the  pseudospectral  algorithm  is  susceptible 
to  apparent  aliasing  instability.  This  problem  , which  was  first  noted 
and  discussed  by  Phillips (1959)  , can  lead  to  erroneous  energy  accum- 
ulation , specially  in  the  short  waves  , and  the  " blow-up  " of  the 
calculation.  We  have  to  remember  that  this  problem  does  not  depend 
on  the  relation  between  the  time  increment  and  gridlength,  and  conse- 
quently is  not  removed  by  simply  decreasing  either  of  these  quantities. 

Merilees(l973  and  1974)  developed  an  algorithm  for  the 
application  of  the  pseudospectral  method  to  the  numerical  integration 
of  the  shallow  water  equations  and  showed  the  ability  of  this  scheme 
to  solve  time  dependent  problems  and  its  superiority  over  the  4th-order 
finite  difference  schemes  to  reproduce  accurate  analytical  solution 
given  the  same  resolution.  Jacques(l976)  used  the  same  methods  to 
integrate  a two-level  model.  In  these  studies  the  problem  of  aliasing 
instability  was  controlled  by  means  of  filtering  where  all  wavelengths 
less  or  equal  to  3-gridlengths  are  eliminated  from  time  to  time. 

Although  using  the  3-gridlength  filter  extends  the  calculation  for  20 
days  without  any  indication  of  instability  , it  is  a crude  closure  appro- 
ximation because  it  is  not  a part  of  the  governing  equations.  On  the 
other  hand  , if  the  Fourier  filter  is  performed  every  time  step  the  model 
would  be  like  a spectral  model  using  Fourier  series  as  basis  functions 
instead  of  spherical  harmonics. 

Jacques(1976)  noticed  that  the  use  of  the  periodic  filter 
produced  substantial  errors  , especially  in  low  resolution  runs.  To 
remove  this  error  and  to  control  aliasing  instability  he  suggested  the 
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use  of  a filter  which  does  not  totally  eliminate  wavelengths  less  than 
three  gridlengths  or  alternatively  the  addition  of  friction  terms. 

In  this  thesis  we  investigate  the  suggestion  of  Jacques  by 
incorporating  a viscous  term  in  the  momentum  equations  and  also 
by  studying  in  more  detail  the  effect  of  periodic  filtering.  A further 
extension  of  the  methodology  results  because  the  viscous  term  requires 
second  derivative  with  respect  to  space  which  previous  applications 
have  not. 

Even  though  Jacques (1976)  was  concerned  with  the  simulation 
of  baroclinic  waves  , the  present  work  is  entirely  concerned  with 
barotropic  motion  as  represented  by  a shallow  layer  of  homogeneous 
incompressible  fluid.  It  will  be  seen  that  there  are  many  interesting 
results  which  accrue  from  these  experiments  which  are  quite  difficult 
to  understand  in  spite  the  barotropic  nature  of  the  model. 


n 
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CHAPTER  H 

DESCRIPTION  OF  THE  MODEL 

II.  1 Basic  equations 

The  equations  of  the  model  are  those  describing  the  flow  of 
a shallow  layer  of  water  on  a rotating  sphere  with  the  addition  of  a 
viscous  force.  One  might  think  of  this  viscous  force  as  an  eddy  diffu- 
sion term  , but  our  purpose  is  simply  to  remove  kinetic  energy  from 
the  system,  especially  from  small  scale  waves  , in  order  to  control 
aliasing  instability.  In  terms  of  spherical  polar  coordinates  ( \ long- 
itude and  latitude  ) the  momentum  equations  are 

/p  VJL  4 .«\_l  . 5 'iVv 


And  the  continuity  equation  is 


• > 


= O 


0) 

w 

U) 


The  advection  and  divergence  terms  are  given  in  spherical  coordinates 

as  - - 


and 


=-Si-  if  + . 

c«i  vf  bA  <*- 

, (Ml  - ' r . \CMc->v)  "I 


where  ^ is  the  velocity  vector  with  component  u and  (east- west 
and  north- south  respectively).  and  Rp  are  the  components  of 

the  viscous  force  with  the  formula  given  in  section  II.  6.1.  The  symbol 
£ takes  the  values  one  or  zero  to  identify  that  the  friction  is  included 


II.  2 Grid  point  equations 

Equations  (1)  to  (3)  are  solved  numerically  on  a latitude- 
longitude  grid  defined  as  follows 

% = 4*0  - T , i ^ T <£  , 

and  = j I <£•  > 

where  - '*Vx  + ^/a.  J = ‘nV^3  a.r\<^ 

4W  and  A*  are  latitudinal  and  longitudinal  increments. 

This  definition  avoids  grid  points  at  the  poles  which  are  considered 
as  singular  points,  as  & \ goes  to  zero.  The  momentum  and  contin- 
uity equations  are  thus  satisfied  every  where  on  the  grid.  If  ^ 
and  ^ denote  the  finite  difference  operators  which  are  to  approx- 
imate X.  , X—  and  -X- 

it  ^ vx 


usual  meteorological  meaning.  With  =0  equations  (l)  to  (3) 

conserve  all  powers  of  potential  vorticity  ( ),  where  ^ is  the 

h 

relative  vorticity.  Also  it  is  possible  to  show  that 

^ \ ^ ^ ) -t-  3 K1,  3 °*fl  = 0 * w 

ft 

where  the  integration  is  carried  all  over  the  earth's  surface.  dA 
represent  an  area  element  ( %A=bX./\Y=st'  cos<f  ).  Equation 

(4)  implies  that 

<f  A 

The  first  term  of  equation  (5)  represents  the  kinetic  energy  while  the 
second  term  represents  the  potential  energy. 


respectively,  then  (l)  to  (3)  can  be 
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written  in  the  form 

= - rav?  ^ 1 Ui 

V— ds*^  «> 

atld  C\\tjW.  + ^SAK  + (,K-a  . ts) 

A centered  finite  difference  scheme  is  used  for  the  left  hand  side  of 
the  above  three  equations,  as  foLlow 

<5  ft  _ fWfc-  ftt-at  a 
* z&t 

where  A stands  for  u,-v  or  h.  For  the  first  time  step,  a forward 
scheme  is  used.  For  space  derivatives  ^ and  , a pseudospe- 
ctral  scheme  is  used  as  described  in  the  following  section. 


II.  3 The  pseudospectral  definition  of  derivatives 

The  definition  of  the  pseudospectral  derivatives  as  used  is 
described  by  Merilees  (1974),  but  for  the  sake  of  completnes s we  shall 
review  the  definition.  In  spherical  polar  coordinates  let  us  consider 
the  variable  A^  defined  at  grid  point  (n,  m)  } A„>w  - A (,An  j'fmV 
*kirt  An=(n-1)&>  ; 1 ^ n ^ 2N, 

and  HJ**  -mt'f  ; 1 ^ m ^ M. 

The  approximation  for  presents  no  special  problem  since  the 

variables  are  periodic  of  A and  do  not  involve  the  pole.  The  lamda- 
derivative  of  An  m is  defined  as 


I v ^-0  iknvX 

^ftk„  = 2T.  , 


mHRi 
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where 


-«kniX 

= -L.  2_  ^n,m  <*,r'd  1 = f~\  * 


*ns\ 


For  the  approximation  of  we  have  to  be  more  careful  since  we 

cross  the  poles.  In  this  case  we  have  to  know  whether  the  variable 
is  vector  component  or  scalar.  A scalar  quantity  is  continuous  at  the 
pole  while,  in  general,  a vector  component  is  not.  The  phi-derivative 
of  is  defined  as 

for  1 fS-  m fS  M and  1 ^ n £*-  N, 


Cw-P 


= n •«*  ^ t'\ 


where 


^ = rn  S“v  ^ * 


For  grid  point  (LM,  LN)  such  that 

LM=N+n  for  1 n N , 

and  LN  = 2M-m+l  for  M+l  m £a.  2M, 

the  phi- derivative  of  A^  w is  defined  as 

= s -£ . -.v  Ul)  e , 


where 


• TR  £ A-.-  e 


In  the  previous  definition  of  the  phi-derivative  s = -l  for  a scalar 

quantity  and  s = +l  for  a vector  quantity. 

II . 4 One  dimensional  semi-implicit  algorithm 

The  numerical  integration  of  the  dynamic  equations  on  a 


j 


latitude-longitude  grids  requires  excessively  small  time  steps  due  to 


F 


r 
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the  convergence  of  meridians,  Merilees  et  al  (1976)  show  that  the  use 
of  space  filtering  in  the  polar  regions  ( generally  used  before  to  avoid 
this  restriction  ) cause  errors  which  lead  to  artificial  transports  of 
mass  and  momentum.  In  the  same  paper  a semi-implicit  technique 
with  one-dimension  is  used  instead  of  the  polar  filter.  The  use  of  a 
semi-implicit  scheme  in  order  to  increase  the  time  step  makes  sense 
from  a physical  point  of  view,  because  it  modifies  only  the  high  freq- 
uency motions  whereas  the  polar  filter  does  not  distinguish  between 
low  and  high  frequency  waves  leading  to  errors  in  the  meteorological 
modes.  The  solution  of  the  one-dimensional  semi-implicit  scheme 
makes  use  of  the  fast  Fourier  transform,  so  the  replacement  of  the 
polar  filter  adds  little  additional  computing  time.  In  equations  (1)  to 
(3)  the  terms  which  give  rise  to  gravity  wave  propagation  along  a lati- 
tude circle  are  treated  semi-in-vplicitly . This  means  that  (2)  will  be 
treated  explicitly.  If  we  define  a variable  A such  that 

= A * ''Pj  •>  ^ 

where  » bX  being  the  time  step,  the  numerical  appro- 

ximation to  the  primitive  equations  (6)  to  (8)  at  grid  point  (I,  J)  and 
time  can  be  written  as 


-oat  T-at 

u.  — m 

2.  at 


nr  ” v 


VT  - 


» ■* 

- ,T*A*  T-  _T_4t 


and 


4M 


KT>‘iKT“*.  j_  r-H*  >"+  ^ kT+  . 


- 1 FT”rt 

4 rv  » 


T+at  ^ T-ft 


Cl®) 

CiO 


Notice  that  and  F^  in  (9)  and  (10)  are  calculated  at  timeT-afc  . 

Equations  (9)  to  (11)  can  be  written  as 


T-f-fct 

U - 


V .!£!  i r*  = v 

* OL  Cos^  A 


- 

CU) 

. ^T,T'kt 

0*> 

. hTJT'a* 

OH) 

T+*t 


The  right  hand  side  of  (12),  (13)  and  (14)  can  be  calculated,  they  are 
function  of  u,  -7  and  h at  times  T*  and  T-At  . This  means  that 

Tfllt  TMt 

can  be  calculated  from  (13).  For  u and  h , let  us  cons- 
ider (12)  and  (14)  for  a particular  zonal  mode  k,  using  the  Fourier 
representation  a variable  A can  be  written  as 


A = 


A Ck)e 


Then  (12)  and  (14)  can  be  written  as 

T+ftt  ,Tt*t 

* W + jrs*  k<K> 


k IK)  V4M 


T,T-at 

= V Ck)  , 

= fiS?  - 


Or 

T+6t 

IMK) 

l 

1 

i 

• 

ik^at 

-I/r,rAt 

1 - vT'T" 

1 

' 

1 

_r" 

» 0*<f  “ il 

..  -T  . 

ikH  at  . 

* c^'f  ' 

which  leads  to 
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V' 
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Mu  ]/[>.»%?]. 

»'■  i {V  " • 


where 


T+*t  X*tK 

Finally  by  inverse  Fourier  transformation  we  get  u and  h 


II . 5 Three  gridlength  and  time  filters 

As  stated  in  the  introduction,  it  is  usually  necessary  to  have 
some  sort  of  smoothing  to  ensure  stability  of  the  calculation  for  more 
than  a few  days.  The  instability  is  believed  to  be  the  so-called  non- 
linear type.  Orszag  (1971)  shows  that  aliasing  can  be  eliminated  if  the 
variables  are  filtered  with  respect  to  three  gridlength  waves  and  sm- 
aller. This  filtering,  used  each  time  step  or  periodically,  for  initial 
conditions  corresponding  to  Haurwitz  waves,  keeps  the  integration 
stable  for  at  least  20  days,  the  maximum  period  of  our  runs.  In  what 
follows  we  call  this  smoothing  the  3-gridlength  filter. 

In  the  3-gridlength  filter,  Fourier  representation  is  used 
for  the  variables  u,-J  and  h then  all  wavelengths  less  than  or  equal 
to  3-gridlengths  are  filtered.  This  filter  is  done  in  both  east-west 
and  north-south  directions.  For  example  with  N points  around  a lati- 
tude circle,  if  An  represents  a variable  before  filtering  and  the 
variable  after  east-west  filtering  then 

j*nA> 

where  . 

*a)  = -V  LAne. 

N nsi 


by  Robert  (1966)  by  the  following  algorithm 

F - F 

T _.T  r T«t  T *T-at 

F = F [F  -ZF  -l-F  J . 

F stands  for  any  of  u,  -1/  and  h.  The  (#)  indicates  a time  filtered 
variable.  In  the  experiments  reported  here  is  set  at  0.02. 


II.  6 Friction  term 

Because  of  the  limited  resolution  of  the  models  which  pre- 


dict geophysical  fluid  motion,  it  is  generally  impossible  to  describe 
the  small  scale  motions.  Through  turbulence  theory,  it  may  be  pos- 
sible to  describe  their  average  behaviour  and  their  effect  on  the  large 
scale  motion.  One  of  the  simplest  of  these  ideas  is  that  of  energy 
cascading  from  large  to  small  scales  to  be  finally  dissipated  by  mol- 
ecular viscosity  at  very  small  scales. 

As  is  clear  from  our  previous  discussion  a friction  term 
will  be  studied  as  an  alternative  to  the  3-gridlength  filter.  The  primary 
purpose  of  this  term  is  to  prevent  non-linear  instability  and,  hopefully, 
to  represent  in  a more  sophisticated  way  what  happens  in  the  atmosph- 
ere. In  the  following  three  sections  a discussion  of  the  viscous  force 


as  included  in  this  work  will  be  presented. 


Mathematical  formula 


In  a general  coordinate  system  the  definition  of  friction  or 
viscous  force  is  given  by  McConnell  (1931)  as 

K - $ > Off) 

O it 

where  0^  is  the  fundamental  or  metric  tensor,  and  3 is  its  contr- 
ol 

avarient  form,  5^,  is  the  Kronecker  delta  and  is  the  viscosity 

stress  tensor.  The  indices  (r,  s,t)  in  the  tensor  notation  take  values 
corresponding  to  the  dimension  of  the  space.  Repeated  indices  are 
to  be  summed  over  the  dimensions  of  the  space  and  a comma  indicates 
covarient  differentiation. 

For  equation  (15)  to  be  useful  we  have  to  express  the  stress 
as  function  of  the  other  variables  of  the  flow.  We  know  that  the  stress 
is  function  of  strain  and  as  usual  we  shall  use  Hooke's  postulate  that 
the  stress  is  a linear  function  of  strain,  or 

= C e-  - 

is  actually  the  velocity  strain  tensor  defined  as 

= "3T  C ^ j 

where  is  the  velocity  vector. 

The  quantities  VrJ  which  form  a mixed  tensor  of  the  forth  order, 
are  called  the  coefficients  of  viscosity.  We  shall  assume  that  the  fluid 
is  homogeneous  and  isotropic.  For  a homogeneous  fluid  the  same 
strain  at  different  point  of  the  medium  produces  the  same  stress.  Thi: 
means  that 
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A necessary  and  sufficient  condition  for  the  above  to  hold  is  that 

wiY\ 

Y =.  o . 

rj,t 

For  an  isotropic  fluid  the  most  general  isotropic  fourth  order 


tensor  with  symmetry  in  (r,  s)  and  (m,  n)  is 


*rs  " / ( ) + A 3 3rs  • 

The  first  term  is  known  as  the  first  viscosity  and  the  second  term  as 

/ 

the  second  viscosity;  and  are  elastic  constants.  It  follows 

that  the  stress  can  be  written  as 

«;,=/*( C ^ + C C ) + A'a'X  , 

= ^ ) + -i:  i 3rs  ( ) , 

= */*  *-rs  + A D > 

= C^s  + ^Sjr  ) + A $rs  ^ ■» 

where  D is  the  divergence.  It  follows  that  (15)  becomes 

Fr»/»a*(Vr,*+Vv)rf  » 

»*[*Wkvr+^+)[a„  , 

where  |<  is  the  curvature  of  the  earth. 

It  is  clear  that  the  last  term  does  not  add  anything  new  to  the  equation 

and  because  we  use  the  viscous  force  to  get  its  smoothing  effect  we 
J 

shall  put  A = 0.  Thus  the  friction  we  have  used  takes  the  form 


While  the  shallow  water  equations  on  a sphere  have  a time-honoured 
usefulness  for  studying  numerical  schemes  and  can  be  derived  with 


. 
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C 


appropriate  limiting  assumptions,  the  form  for  a spherical  friction 
force  has  not  been  so  carefully  treated  especially  as  to  its  limiting 
form  for  quasi-horizontal  flow.  We  have  arbitrarily  chosen  to  consi- 
der the  flow  to  be  strictly  two  dimensional.  However  it  will  be  our 
contention  that  a simple  friction  term,  with  coefficient  of  viscosity 
independent  of  the  flow,  is  not  very  useful  for  controlling  aliasing  ins- 
tability for  the  simple  reason  that  the  Reynolds  number  required  at 
the  scale  of  the  grid  is  simply  too  small  not  to  affect  the  scales  of  int- 
erest. In  other  words  the  magnitude  of  friction  required  to  control 
numerical  instability  excessively  smooths  the  scales  of  interest. 

For  surface  spherical  polar  coordinate  in  two  dimension 
> and  tp  ; tp  is  the  latitude  and  ^ is  the  longitude.  The  distance 
ds  is  defined  as 

(cJs)  - cl'CcJ'O  + cc  (o/,A)  J 

= + \x  • 

So  ctnd  ^ c*'x\ 

3"  = i/o^  , 

where  a is  the  radius  of  the  earth. 


As  shown  in  appendix  (I)  we  find  that  and  are  given  by 

F - $ 5 i±.  +-i-  _3S!^  4 » W 


— ~ ^ O-  -W  — — ^7 

v ^ J + \X' 


ecnJl 


BHMMMHMP 
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, J_  W.  . _1_  irL  +3WJ  w 
Co\y  v*U  ^ o.»'y  *XV  \J 

T ^ c*i  1 


In) 


where  u and  -0  are  the  horizontal  components  of  the  wind  vector  and 
$ is  the  kinematic  viscosity  coefficient. 


II. 6. 2 Analytical  effect  of  friction 

One  of  the  useful  features  of  the  friction  term  used  is  that 
if  the  flow  is  ( horizontally  ) non-divergent,  then  it  can  be  shown  that 

IK  is  the  unit  vertical  vector,  'f  is  the  horizontal  friction  vector 
and  J is  the  relative  vorticity.  Since  the  condition  of  non-divergence 
implies  that  the  shallow  water  equations  reduce  to  the  vorticity  equation, 
then  using  the  previous  formula  and  considering  that  the  change  in 
vorticity  is  due  to  friction  we  can  write 


Now  if  we  write  the  vorticity  as  a pure  spherical  harmonic 

, ImX 

l = A<t)  e.  , 


where  A(t)  is  the  amplitude  of  the  vorticity  which  is  function  of  time, 

£ (*)  is  the  associated  polynomial  of  the  first  kind,  m is  the  east-west 
wavenumber  and  n is  the  two  dimensional  wavenumber,  then 


fwf*  if 


-£ "]  Atfl  P„«>  £ * 


I 


r 
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or 


AU) 


= ~lrln0'+,)"2'^  • 


Integrating  we  get  A ■=•  A0  » 
where  *<•  = ^ ^Cn+0-2.^  . 

So  we  can  see  that  the  effect  of  friction,  is  to  reduce  the  amplitude  of 
the  vorticity  exponentially  with  time.  The  effect  will  be  larger  for 
small  wavelengths  or  large  wavenumbers,  as  it  is  clear  from  the  for- 
mula for  . 


II.  6. 3 The  effect  of  the  friction  term  on  divergence 

The  initial  conditions  used  in  our  experiments  are  such 
that  the  divergence  and  its  first  time  derivative  is  initially  zero.  How- 
ever, due  to  non-linear  interaction,  even  with  the  above  initial  condi- 
tion significant  divergence  develops.  The  problem  of  aliasing  is 
associated  with  the  appearance  of  wiggles  and  increase  of  divergence. 

In  this  section  it  is  desirable  to  show  that  the  amplitude  of  the  diverg- 
ence field  will  decrease  due  to  the  existence  of  the  friction  term,  or 
in  other  words  the  divergence  field  is  as  well  controlled  by  the  friction 
force . 

Using  equations  (16)  and  (17)  with  the  definition  of  divergence 
in  spherical  polar  coordinate,  it  is  possible  to  show  that 


v.  F = 


\ 

ot  OxS'f 


[ 


*L  Y_(K,o*Q 

IX  v? 


J 


F 
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or 


v*  , ■»  3 . ^ 

TJ . |F  = vj_  _i 4.  *-  - kL  v + _£ \it_ 

a?  'bvs’  a*  c.^  ^'iX  ^ oc*  c.*4*  VOX  o?C.i\  *X* 


1 "X*,  _ 4.  V^*f  Eli.  (i+Jl\*± 

^»a ?V  *OX  T <*?c.»*«  'b'f’-  iu'ew'f  v>f  Vc-r 


a? c.!v<f  v**>  ’ <*-  C*»*$  ’b'T  b>4  oc'c-f 

^ \-J  2.  /-•_.*.  ****L 

•?C-a'Vf  Vf 


*x 
00 


and 


rj*Vv  _ \ . J ^ , _\__  j^L  __  a***  ibL 

ou*C>  *>v  «?GuVf  o>  **"*  “Vc*1*  VXII. * * V 


q.  Sw*f  \v*~  1 \ \ Ub  ■any*' 

"^“  oi  Cotv  Vf\>  o?<u*'i  "Vf^X  o?c»»'f 

-l  / M^-y\  H ._  v -J  /.  «\ 

So  from  (18)  and  (19)  one  can  conclude  that 

^.\F=  a.C'Q+Jf)  fc 

If  we  take  the  divergence  of  the  momentum  equation  and  considering 
that  the  local  change  of  the  wind  is  governed  only  by  the  friction  we 
get  the  formula 

J£  = i*(vVV)^ 

From  the  above  equation,  it  is  clear  that  the  friction  reduces  the 
amplitude  of  the  divergence  except  probably  for  the  very  large  scales. 


S'^tf  V~-J 


C 


II.  7 The  criteria  for  numerical  stability 

The  guidelines  for  the  choice  of  time  step  for  the  integration 
of  equations  (1)  to  (3)  have  been  determined  by  an  analysis  of  the  CFL. 
criteria  in  Cartesian  geometry  as  given  in  appendix  II.  These  criteria 


are  summarized  below. 

1.  Neglecting  the  friction  terms  and  using  a fully  explicit  scheme, 
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C 


the  condition  is 

at^  l/<LvioV 

— H. 

where  c = (gH)  is  the  speed  of  gravity  waves;  L and  K are  respect- 
ively the  maximum  wavenumbers  permitted  in  the  north-south  and 
east-west  directions.  As  shown  in  appendix  II  the  inclusion  of  friction 
terms  in  fully  explicit  scheme  does  not  affect  the  CFL.  criterion 
significantly  in  our  case. 

2.  Neglecting  the  friction  terms  and  using  the  one  dimensional 
semi-implicit  scheme  as  described  in  n.4  the  condition  reads 

At  i/l—C.  (2.1) 

Due  to  the  convergence  of  meridians  near  the  poles  K^L  and  thus 
the  one  dimensional  semi-implicit  algorithm  permit  much  longer  time 
step  than  those  given  by  the  fully  explicit  scheme.  However,  if  two 
dimensional  semi-implicit  scheme  is  used  all  the  modes  will  be  neutral. 

3.  With  friction  and  the  one  dimensional  semi-implicit  algorithm 
we  have  two  conditions  to  be  satisfied,  they  are 

at  ^ 1 A>> ( kV C)  • i«> 

i & $ Wx £ + * ut)  (kVC),  (j.%) 

=gH. 

In  fact,  we  can  summarize  the  result  of  this  case  as 

a.  If  $ is  relatively  small  the  condition  for  stability  is  (21). 

b.  As  ^ becomes  larger  the  condition  for  stability  goes  to  (22). 
To  show  the  difference  between  the  time  steps  as  given  by  the  conditions 


and 

where 


Experiments  with  3-gridlength  filter  applied  each  time  step  permit 
a longer  time  step  due  to  the  removal  of  short  waves. 


II. 8 Computing  time  per  step 

A CDC  7600  computer  in  single  precision  arithmetic  ( 60 
bit  word  length  ) was  used  for  all  the  experiments  discussed  in  this 
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thesis.  The  time,  in  seconds,  required  per  step  of  integration  for 
the  one  level  model  used  are  approximately  as  follows: 

a_  Including  Fourier  filter,  applied  to  u,  V and  h,  for  two 
successive  steps  each  three  hours 
T = 0.22x10*  D. 

b.  Including  Fourier  filter,  applied  to  u,  V and  h,  each  time 

step 

T = 0.27*103  D. 

c.  Including  the  viscous  force 

T = 0.32*10*  D. 

d.  Including  neither  the  viscous  force  nor  the  Fourier  filter 

T = 0 . 21  x 10*  D. 

D is  the  number  of  degrees  of  freedom  of  the  model. 

This  means  that  the  use  of  the  periodic  filter  ( two  success- 
ive steps  per  three  hours  ) increases  the  time  of  calculation  by  less 
than  5%  ; filtering  each  time  step  by  less  than  29%  ; while  the  viscous 
force  increases  the  time  by  about  52%. 


■ »«- 
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CHAPTER  IH 
INITIAL  CONDITIONS 


III . 1 Initial  condition  for  the  test  experiment 

The  field  is  initially  non  divergent  or  u and  -V  are 
calculated  from  a stream  function  defined  as 

= <£  K (£•*  f)  C.i  C"'A”wt ) - « K J * ot) 


The  height  field  is  given  by 

K S ^ 0*1  'f)  (Sn^'  wut)  J Cj-S") 

3 3 

where 

m is  the  wavenumber, 
t time. 


k,  k,_  , ft  , w and  b^are  constant, 

k,  = k eT*  , OA) 


•<  is  function  of  the  wavenumber  ( see  11.6.2  ), 

-ft-  is  the  angular  velocity  of  the  earth, 
u,  -V  and  h given  by  (24)  and  (25)  do  not  generally  satisfy  (1)  to  (3).  To 
make  u,  v and  h satisfy  (1)  to  (3)  for  all  times,  an  extra  term  is  added 
to  each  of  the  equations  (1)  to  (3).  The  modified  equations  can  be 
written  as 


and 


u 

f — 

T <*■  Cm'S 

*X 

Vi 

_l.  _ih_ 

VI 

>t 

^OL  Cal'? 

1L 

_r 

u 

a Cos'f  L. 

— v~r  — £ nj  -v  -5 — - FT_  G*  = o , 


Tx 


*X  ^ J 


CM) 
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C^) 


* i.e.  the  wave  component  is  proportional  to  » A'?) 


- — ' ■ • " 
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where  (*  = t + u tan  vf  /a. 

The  expressions  for  Gu,  Gv  and  Gh  are  obtained  by  substituting  u.  -V 
and  h from  (24)  and  (25)  in  (27)  and  (29).  This  gives 

6Jul  r (tm  CK(- Kt)  + •». 

— «ck»  C»sC«>-ut) ^ -«■'*' k,(c»nO  £ l4C»+"0  - 

4-o-  Sif\^^X-u^)  + a<-  CoS  A-*-i ^ ^ 

— 2_  *a(iavO  k*(S.o*4)  <oJt**A-at)  £c***0  c-*^ -**3  * 

04  =■  SiA^(c»s^)  U k,  C.if.tnA-*-1^)  + <*«<  »*  k»  SinC'+A-'-tk) 

_ a.  w\l  f\  k,  to»W-^0  + a.  JV.  kv«-  *a  C*sC*^A-*Jt) 

4-  2a.  A kvrr\  Cos  ^ t)  — Z.  -d_c«.  kv  ^ CttsC*'A-‘Jfc)3 

4-S»V4  (Coi4)W'+l  kxo>*»\  G>i^A-wt)_  20. a k4  C-s^A-wt) 

— 2.0. a )<\  tosC»*w>-ut)  4-  Coj c^A-^*) 4-  aJL  ok^ Cosd^A-^f )J 

(Sojy)  3 [~a-*^  k,1" _». «\  kt6*-0  SuvC^A- 

4.  <*.  wv-  k'  CosV (.*->-«-» t)  4-  ot  kn  kl-  s« A C*aA -*»*)+•  (*"-0  k,  $«  * 6"*>-u  0 

— 20L  k^K1  Cos\*aA-«0  _ ia  C»A  Cos’CwsA-ut)]] 

+ s»ny  ^jT* £"•  ^ * c**0*A -°t)  + k^cJ<l*AA-«t) 

4.  a.  k,  C«s  4-  a.  A ?«kV{  C.»v^ 

_2Ji*  mx(*vt-l)  5*a4  (Co*  *?)'*’  SiA(.^A-ot)  > 

and 

£k  = ^r^°L-  ^in^A-ut)  jivu+A»  (k,-kv)-»Akt k,C-,r 

x c^vf)1  cco^r. 


r 


With  no  friction  ( ^ =0)  the  solution  of  the  system  (27)  to  (29),  with  the 
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initial  condition  given  by  (24)  and  (25),  represents  a wave  which  moves 
with  angular  phase  speed  u>  without  change  in  shape  or  amplitude. 
However,  with  friction,  the  wave  will  move  with  the  same  phase 
speed  without  change  in  shape  but  the  amplitude  for  u and  V will 
decrease  according  to  the  relation  A = A,  g.  , where  A is  the 
amplitude  at  time  t and  A,  is  the  initial  amplitude.  In  fact,  this  is 
why  we  put  the  relation  (26),  because  we  know  from  section  II.  6. 2 
the  effect  of  friction.  The  values  used  in  the  test  experiment  are 
m = 4 , ( wavenumber  4), 
k = = A =7 . 848 x 10  ^ sec', 

-rt.  = 2TT/86400  sec', 
h#=  3 km, 
a = 6 . 4x 10*  km, 

3 = 10  m/sec*. 

HI.  2 Haurwitz  wave  initial  condition 

The  initial  field  is  non  divergent,  so  that  u and  v can  be 
derived  from  a stream  function.  In  fact,  we  use  the  same  stream 
function  of  the  test  experiment,  or 

set  k C.C.^'T  ***  . 

The  height  field  which  balances  this  initial  flow,  or  the  height  field 
which  makes  the  time  derivative  of  divergence  initially  zero  is  given 
by 

ft  Of)  + C.I  RX  4 * &(*) C.* * 


I 


I 


where 


V--  - - 
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AW)  = "3 

M0-  c~[<„V"'^)-^+'>t‘'l. 


bC'f'i  _ _K_l  <y*  £(wn+i)  cl  - On-*-*)"!  , 


and  C = cos(tf)  . 

The  values  used  in  the  experiments  are 

m = 4 and  6,  ( wavenumbers  4 and  6 ), 

•t  —l 

k = f\  = 7.848x10  sec*  , 

-A-=  2TT /86400  sec"1  , 
h,=  8 km  , 
a = 6 . 4x  10*  km  , 

3=10  m/sec1 . 

Although  the  divergence  and  its  time  derivative  is  initially  zero, 
significant  divergence  develops  and  thus  the  waves  do  not  propagate 
with  the  Rossby-Haurwitz  phase  speed,  nor  do  they  mantain  their 
shape  precisely. 


r 

j 
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CHAPTER  IV 

EXPERIMENTS  AND  RESULTS 


c 


A basic  mean  flow,  independent  of  longitude,  superposed 
on  it  a disturbance  will  be  integrated  with  time  using  equations  (1)  to 
(3).  Two  different  types  of  disturbance  are  used,  wavenumber  four 
whose  form  changes  a little  during  more  than  fifteen  days  of  integr- 
ation (stable  wave),  and  wavenumber  six  which  breaks  down  completely 
within  few  days  forming  cut-off  lows  and  a region  of  zonal  easterlies 
at  about  50*  N (unstable  wave).  Experiments  with  the  following  diff- 
erence will  be  described 

a.  Runs  with  Fourier  filtering. 

b.  Runs  including  the  viscous  terms  in  the  momentum  equations. 

c.  Runs  which  do  not  include  either  of  the  damping  forces  (a) 
or  (b). 

IV.l  Experiment  notation 

Hereafter  we  shall  use  the  experiment  notation  AHJJ(L-M), 

where 

A stand  for; 

F if  friction  is  used, 

P if  Fourier  filter  is  used, 

N if  neither  friction  nor  Fourier  filter  is  used, 

H is  the  number  of  grid  points  per  latitude  circle, 

JJ  is  the  number  of  grid  points  between  the  poles, 


HM 
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L will  be  used  in  case  of  friction  and  Fourier  filter.  In  case  of 

friction  it  indicates  the  value  of  kinematic  viscosity  used  ($  ),  e.g. 

. Hr  1. 

experiment  including  friction  with  v = 10  m /sec,  L = 4.  Experiments 
with  Fourier  filter  have  L to  indicate  the  value  of  the  number  of  succ- 
essive time  steps  used  per  M hours.  If  L and  M are  not  included 
in  an  experiment  notation  for  Fourier  filter,  this  means  that  the 
3-gridlength  filter  is  used  each  time  step. 

For  example  in  an  experiment  with  11  = 32  and  JJ=16, 

the  experiment  notation  is  F3216(5)  if  friction  is  used  with  0=10  rn/sec, 
the  experiment  notation  is  P3216  if  the  Fourier  filter  is  used  each  step, 
the  experiment  notation  is  P3216(2-3)  if  the  Fourier  filter  is  used 
for  two  successive  steps  each  three  hours, 
the  experiment  notation  is  N3216  if  neither  friction  nor  Fourier 
filter  is  included. 


In  (A 


30 

TIME  (hours) 


The  logarithm  of  the  square  of  the  amplitude  of 
Fourier  wavenumber  4 of  the  V- component  of 
wind  (summed  for  all  grid  latitudes)  as  function 
of  time  for  the  test  experiment,  F3216(7). 
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IV.  2 Test  experiment 

It  is  essential  before  starting  our  experiments  to  be  sure 
that  the  effect  of  including  the  viscous  force,  in  the  momentum  equations, 
gives  the  expected  results.  In  other  words,  it  is  important  to  test  the 
program  with  and  without  this  force  for  a simple  initial  condition  to  be 
sure  that  there  are  no  programming  errors  or  hidden  instabilities. 

The  initial  condition  used  in  this  test  is  Haurwitz  initial 
condition  for  wavenumber  four  with  the  forcing  function  included  as 
discussed  in  section  III.l.  This  means  that  the  model  represent  a 
non-diver-'ent  barotropic  fluid  for  all  times  of  integration.  This  means; 

a.  If  the  model  runs  without  the  viscous  force,  the  result  must 
be  just  a moving  wave  with  known  phase  speed  and  without  change  of 
amplitude . 

b.  If  the  model  includes  the  viscous  force,  then  the  u and  V 
fields  will  move  with  the  same  phase  speed  but  their  amplitudes  will 

decrease  exponentially  with  time  according  to  the  equation,  (see  II. 6. 2) 

-•tt  v 

A = , where  «*<=')  m(m-r3)/a. 

The  previous  equation  can  be  written  as 

log  (A)X  = log  (Aa)X - 2<xt.  (30) 

This  means  that  the  logarithmic  change  of  the  amplitude  squared  with 
time  is  linear  with  a slope  of  -2o<.  In  this  test  experiment  \)  was 
taken  to  be  l(?  rn“ /sec  giving  an  e-folding  time  of  roughly  three  days. 
The  result  for  the  experiment  F3216(7)  is  shown  in  Fig.l.  As  expected 
from  formula  (30)  a straight  line  resulted  on  a graph  with  the  coordinate 
log  (A)  and  t.  From  the  graph  the  slope  of  the  line  2 =0.137x10 
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and  from  the  values  used  2 ®<  =0.1376x10  . In  this  experiment  the 
time  step  used  is  2 minutes.  There  was  no  evidence  of  any  numer- 
ical difficulties. 


IV.  3 Experiments  with  wavenumber  four 

Hoskins  (1973)  has  shown  that  a Haurwitz  wavenumber  four 
superimposed  on  a constant  angular  velocity  zonal  flow  is  stable  in 
the  sense  that  the  wave  propagates  in  an  east-west  direction  with  little 
change  in  shape.  In  other  words,  there  is  little  cascade  of  energy  or 
enstrophy  during  a time  integration.  As  such,  one  expects  that  there 
will  not  be  a serious  difficulty  in  the  simulation  of  this  wave  from  the 
point  of  view  of  the  control  of  aliasing  instability.  Indeed  we  have 
found  that  an  integration  with  low  resolution  (11  = 32  and  JJ=16)  can  pro- 
ceed for  16  days  without  any  controls  at  all.  On  the  other  hand,  the 
inclusion  of  filtering  or  a friction  term  while  permitting  a longer 
integration  time  has,  surprisingly,  rather  significant  effects  on  the 
nature  of  the  solution.  In  this  section  we  describe  the  result  of  expe- 
riments with  low  resolution  and  high  resolution  (11  = 64  and  JJ=32).  We 
found  that  the  main  wave  moves  eastward  with  faster  speed. for  high 

resolution  experiments  . Its  displacement  varied  between  10.7*per 

0 o 

day  and  11.4  per  day  compared  to  the  12.2  per  day  predicted  by  the 
nondivergent  theory.  The  phase  speeds  obtained  here  do  not  contradict 


* As  shown  by  Haurwitz  (1940),  in  a nondivergent  barotropic  atmos- 
phere the  flow  pattern  moves  from  west  to  east  without  change  of  shape 
with  the  angular  velocity  v,  where 

v=[R(R  + 3)W  -2-nJJ/(l-rR)(2+R),  the  symbols  as  in  HI.  2. 
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LONGITUDE 


Initial  conditions  (Haurwitz  wave  4);  u(m/sec) 
v(m/sec)  and  geopctential  height  h(meters) 
with  interval  240  meters. 
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those  obtained  by  Bourke(1972)  and  Merilees  (1974),  as  some  of  the 

constants  are  different  with  the  result  that  in  these  experiments  we  get 
faster  waves.  The  phase  speeds  obtained  here  are  very  similar  to  the 
values  obtained  by  Hoskins  (1973),  while  the  values  obtained  by  Phillips 
(1959)  are  smaller,  perhaps  due  to  the  use  of  second  order  finite 
differences.  The  initial  conditions  of  u,  v and  h are  shown  in  Fig.  2. 

IV.  3.1  Low  resolution  experiments 

We  compare  the  results  of  experiments  F32l6(5),  P32l6(l-3) 
andN32l6.  With  this  resolution  the  basic  disturbance  is  described  by 
8 grid  points  along  a latitude  circle.  The  maximum  zonal  wavenumber 
which  is  completely  resolved  is  15.  The  filter  used  in  the  experiment 
P32l6(l-3)  will  periodically  remove  all  wavenumbers  greater  than  10. 
The  time  step  used  in  these  experiments  is  10  minutes. 

In  each  of  the  above  mentioned  experiments,  the  behaviour 
was  substantially  the  same  for  the  first  3 or  4 days.  The  disturbance 
quickly  developed  a north-west  to  south-east  tilt  in  the  northern 
latitudes  of  the  northern  hemisphere  and  a north-east  to  south-west 
tilt  in  the  polar  latitudes.  This  has  the  effect  of  transporting  westerly 
angular  momentum  from  the  northern  latitudes  towards  both  the  equator 
and  the  pole  (see  appendix  HI).  At  this  point,  the  behaviour  of  the 
disturbance  as  indicated  by  experiment  P3216(l-3)  begins  to  differ 
from  the  other  two  experiments.  In  the  experiment  P32l6(l-3)  the 
wave  tilt  developed  initially  is  maintained  whereas  in  the  other  exper- 
iments the  wave  tilt  reverses.  These  results  are  shown  in  Fig.  3(A) 
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Latitudinal  distributions  cf  phase  angle  (A)  and  lon- 
gitudinal zonal  velocity  (3)  on  the  days  3,  o ^r-d  9 
(from  left  to  right)  for  the  experiments  ?32l6(l-3) 
(solid  line),  ?3 216(5)  (dash  line)  and  N3216  (dash- 
dot  line).  Dash-cross  line  used  when  the  experiments 
ia  fnr  initial  nrof  11 
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where  we  have  plotted  the  phase  of  the  main  disturbance  (wavenumber  4) 
for  the  V-field  as  a function  of  latitude.  The  figure  indicates  that  the 
behaviour  of  the  flow  as  described  by  experiments  N3216  and  F3216{5) 
is  that  of  a stable  nature,  that  is,  tending  to  return  to  the  initial 
condition.  On  the  other  hand  in  experiment  P3216(l-3),  the  zonally 
averaged  u-component  continues  to  accelerate  in  low  latitudes  and 
decelerate  in  middle  latitudes.  As  a result  the  distribution  of  the 
zonally  averaged  u-component  of  the  wind  is  completely  different  after 
9 days  of  integration  as  indicated  in  Fig.3B.  As  we  shall  show 
presently,  the  behaviour  of  the  experiment  P 3 216(1-  3 ) appears  to  be 
anomolous  and  due  to  the  filtering  applied. 

IV  . 3 . 2 High  resolution  experiments 

Since  we  expect  that  the  higher  resolution  experiments 
should  provide  a more  accurate  approximation  to  the  true  solution,  we 
are  able  to  use  such  experiments  to  tell  which  of  the  low  resolution 
experiments  are  going  wrong.  In  this  section,  we  compare  experiments 
N6432,  F6432(5),  P6432(2-3)  and  P6432  for  the  same  initial  conditions 
as  in  the  low  resolution  experiments.  In  this  grid  system  the  maximum 
zonal  wavenumber  which  can  be  completely  represented  by  the  grid  is 
31.  The  filter  when  applied,  will  remove  wavenumbers  greater  than 
20.  The  time  step  used  in  these  experiments  is  four  minutes. 

* Experiment  P6432(2-3)  is  discussed  and  not  P6432(l-3)  because, 
as  we  shall  see  in  section  IV.  3 . 3 experiment  P6432(l- 3)  exites  the 
two  time  step  computation  mode. 
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Fig.  4.  Same  as  Fig.  3(A)  except  for  high  resolution 
experiments . 
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Fig.  5,  Isopleths  for  the  u-field,  oc  day  9,  for  the 

experiments  F3216(5).  N3216,  P32l6(2-3),  F6432(5) 

and  P6432(2-3);  A,  B,  C,  D and  E respective! 
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Fig.  4 shows  the  distribution  of  the  phase  angle,  for  the 

i 

V-field,  against  latitude  for  the  experiments  F6432(5),  N6432  and 

P6432(2-3)  after  3,  6 and  9 days.  Fig.  5(A,  B,  C,  D and  E)  shows  the 
u-field  on  day  9 for  the  experiments  F3216(5),  N3216,  P3216(2-3), 

F6432(5)  and  P6432(2-3).  It  is  possible  to  conclude  from  the  results 

t 

of  high  resolution  experiments  compared  to  low  resolution  experiments 
that; 


1.  The  high  and  low  resolution  models  are  in  good  agreement  up 
to  3 days.  In  fact,  the  fields  of  u,  V and  h are  about  the  same  up  to 

3 or  4 days  for  all  the  experiments  and  thus  they  are  not  shown. 

2.  The  high  resolution  models  generally  give  the  same  results 

up  to  9 days.  There  are  some  differences,  the  most  noticible  of  which 
is  shown  in  Fig.5(D  and  E)  where  the  low  centre  in  the  equatorial  area 
for  the  u-field  has  minimum  value  of  westerly  wind  of  the  order  of 
10  m/sec  for  F6432(5),  while  the  same  centre  has  a value  of  about 
30  m/sec  for  P6432(2-3).  The  difference  for  the  u-field  for  low 
resolution  experiments  after  9 days  is  more  pronounced  as  is  clear 
from  Fig.  5 (A,  B and  C). 

3.  In  low  latitudes  the  wave  moves  faster  in  the  high  resolution 
experiments,  further  there  is  less  difference  in  phase  speed  between 
low  and  high  latitudes  in  the  case  of  high  resolution  indicating  that  the 
wave  is  more  stable. 

4.  The  change  of  phase  tilt  with  time  in  high  resolution  experiments 


I 


i 


**  For  the  same  reason  as  the  footnote  of  page(34),  experiment 
P32l6(2-3)  is  discussed  and  not  P32l6(l-3). 
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is  indicative  of  a stable  flow.  Fig.  4 shows  that  the  momentum  trans- 
port on  the  third  day  is  southward,  on  day  6 is  northward  from 
low  latitudes  while  there  is  a weak  southward  transport  of  momentum 
on  day  9.  On  the  other  hand  in  the  low  resolution  experiments  with 
the  periodic  filter,  momentum  is  continuously  transported  from  middle 
latitudes  to  the  north  and  south. 

Fig.6(A  and  B)  shows  the  distribution  of  the  mean  flow  for 
the  u-field  against  latitude,  on  day  9,  for  the  experiments  F6432(5), 
P6432(2-3)  and  P6432.  F6432(5)  and  P6432  give  the  same  distribu- 
tion but  P6432(2-3)  gives  a weaker  mean  flow  for  the  middle  latitudes 
and  higher  values  for  the  north  and  south.  It  is  possible  to  conclude 
from  what  we  discussed  that  the  behaviour  of  the  flow  for  P32l6(2-3) 
is  incorrect  due  to  the  application  of  the  periodic  filter.  On  the  other 
hand,  the  high  resolution  experiment  with  periodic  filtering  is  able  to 
capture  the  oscillation  in  momentum  transport  which  leads  to  a stable 
wave.  We  conclude  that  some  of  the  waves  that  are  filtered  in  the  low 
resolution  experiments  are  essential  components  of  a stable  oscillation. 

When  no  smoothing  is  applied  to  the  high  resolution  exper- 
iments the  calculation  " blows-up  " after  about  8 days.  Note  the  noise 
on  the  distribution  of  the  mean  flow  as  shown  in  Fig.  6(A).  We  can  also 
see  irregularities  in  low  latitude  in  the  u-field  as  shown  in  Fig.  5(D). 
The  problem  of  instability  will  be  discussed  further  in  section  IV. 3. 4. 


m/sec 
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TIME  (hours) 


TIME  (hours) 

The  amplitude  of  wavenumber,  12  of  the  "7- component  of 

wind  as  a function  of  time  at  two  latitudes.  Solid 

lines  are  for  5.6*N  and  dotted  lines  for  39.3*  N.  The 
upper  figure  (A)  corresponds  to  the  experiment  N3216 

the  middle  (B)  to  P3 216 (1-3)  while  the  lower  (C)  to 
P32l6(2-3).  The  circles  on  (b)  represent  the  experi- 
ment without  Robert  time  filter. 
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IV.  3. 3 Further  experiments  with  low  resolution 

Because  of  the  anomalous  behaviour  of  the  low  resolution 

experiments  with  the  filter  it  was  decided  to  investigate  in  greater 
detail  the  effect  of  the  filter.  The  first  step  in  this  process  was  to 
study  the  variation  in  time  of  one  of  the  wavenumbers  before  and  after 
a filtering  time.  We  choose  to  look  at  the  behaviour  during  the  first 
few  hours  of  the  integration  both  because  it  is  more  convenient  and  it 
should  indicate  the  tendency  for  the  wavenumber  four  flow  to  cascade 
energy  and  enstrophy.  Fig.  7(A)  illustrates  the  behaviour  of  the  amp- 
litude of  wavenumber  12  in  the  V-field  when  no  filtering  is  applied 
while  Fig.  7(B)  indicates  what  happens  when  filtering  is  applied  at  3 
hours  from  initial  time.  We  have  plotted  these  curves  for  two  latitudes, 
one  near  the  equator  and  the  other  in  middle  latitudes,  we  can  see 
that  the  behaviour  is  quite  similar.  We  note  that  the  application  of 
the  filter  at  one  time  step  has  the  effect  of  exciting  a two  time  step 
computational  mode.  Subsequently  this  mode  is  damped  out,  but  much 
faster  than  would  expected  from  the  effect  of  the  Robert  time  filter, 

Fig.  7(B).  It  is  quite  simple  to  understand  why  this  occurs  especially 
if  the  time  derivative  of  the  amplitude  of  this  wave  is  essentially  inde- 
pendent of  its  own  amplitude.  In  fact,  the  subsequent  behaviour  of  the 
wave  amplitude  is  consistent  with  the  idea  that  there  is  some  quasi- 
steady amplitude  which  should  be  acheived  in  this  wavenumber.  If  that 
were  the  case,  we  should  expect  to  see  the  amplitude  to  increase  when 
it  is  excessively  low  and  remain  about  the  same  level  where  it  is  before 
filtering.  These  ideas  are  further  supported  by  the  curves  in  Fig.  7(C) 
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which  show  the  result  of  applying  the  filter  on  two  successive  time 
steps.  Note  how  the  two  time  step  oscillation  is  considerably  reduced. 
Further,  the  amplitude  of  the  wave  tends  to  grow  back  to  its  level 
before  filtering.  This  indicates  that  energy  will  tend  to  flow  to  this 
wave  when  it  has  little.  Since  it  is,  in  general,  good  practice  to  mini- 
mize the  exitation  of  two  time  step  oscillation  it  was  decided  that  it 
would  be  better  to  apply  the  filtering  on  two  successive  time  steps. 

One  of  the  rather  arbitrary  features  of  the  application  of 
the  filter  is  how  frequently  it  is  applied.  Some  very  rough  experiments 
by  Merilees  (1974)  suggested  that  every  three  hours  was  alright,  but 
no  serious  effort  was  made  to  investigate  it  more  thoroughly.  Here 
we  carried  out  a number  of  experiments  to  see  how  the  solution  changes 
as  we  change  the  frequency  of  the  application  of  the  filter. 

We  first  compared  the  effects  of  filtering  at  two  successive 
time  steps  every  three  hours  and  every  hour  with  filtering  once  every 
three  hours.  Fig.  8(A)  shows  the  distribution  of  the  phase  angle  of 
zonal  wavenumber  four  of  the  V-field.  We  note  that  the  tilting  of  the 
phase  is  less  pronounced  when  we  filter  every  hour.  Filtering  at  two 
successive  time  steps  every  three  hours  tends  to  produce  slightly  less 
tilt  but  the  distribution  of  phase  is  quite  similar  to  that  obtained  by 
filtering  once  every  three  hours.  In  Fig.  8(B)  we  show  the  zonal  com- 
ponent of  the  u-field  after  9 days  and  we  can  see  the  effects  of  the 
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* This  practice  is  apparently  well  known  by  those  who  have  experience 
in  grid  point  models.  However,  people  who  work  with  spectral  models 
have  generally  not  concerned  themselves  with  it. 


PHASE  ANGLE  Idagreea) 


A:  The  phase  angle  of  wavenumber  4,  for  the  v-field, 

as  function  of  latitude  after  3,  6 and  9 days  (from, 
left  to  right).  B:  The  corresponding  zonal  average  of 
the  u-field  after  9 days. 
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Fig.  9(A  and  B).  Same  as  Fig.  S(A  and  B) 
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Fig.  10.  Same  as  Fig. % but  for  the  experiment  P3216. 
Dotted  line  is  the  initial  profile. 
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tilting  of  the  phase  lines  . While  there  are  some  differences  caused 
by  the  double  application  of  the  filter  at  three  hour  intervals,  the 
interval  of  time  between  the  application  of  the  filter  appears  to  be  more 
important.  Based  on  this  consideration  a number  of  experiments 
were  performed  in  which  the  time  interval  between  application  of 
filter  was  varied  between  zero  and. 48  hours.  In  each  of  these  experi- 
ments the  filter  was  applied  at  two  successive  time  steps.  Interest- 
ingly, we  find  that  the  tilting  of  the  phase  of  the  wavenumber  four  of 
the  V-field  is  not  monotonically  dependant  on  the  time  interval  between 
application  of  the  filter.  In  fact,  the  application  of  the  filter  at  12 
hours  intervals  appears  to  produce  the  maximum  effect  as  shown  in 
Fig.  9(A).  Naturally,  a similar  effect  is  produced  on  the  profiles  of 
the  zonally  averaged  component  of  the  u-field  as  shown  in  Fig.  9(B). 
Using  the  filter  each  time  step,  surprisingly,  showed  the  more  stable 
solution,  Fig . 10(A  and  B),  with  the  least  change  of  energies.  Fig. 11 
shows  the  kinetic  and  potential  energies  as  well  as  their  sum  during 
a period  of  13  days  for  the  low  resolution  experiments.  The  figure 
shows  a change  of  less  than  2%  in  the  kinetic  energy  for  the  experiment 
P3216  during  the  13  days  of  integration. 

As  a convenient  summary  of  the  sort  of  results  obtained 
we  present  Fig.  12.  The  ordinate  is  the  difference  between  the  zonally 
averaged  u-field  after  9 days  of  integration  and  the  initial  value  of 
this  same  field  at  about  51* N.  Since  the  true  solution  oscillates  about 
the  initial  conditions,  this  provides  a measure  of  error.  Along  the 
abscissa  we  plot  the  different  experiments  in  order  of  decreasing  time 


Fig.  li.  The  percentage  change  from  initial,  time  of  the 

Kinetic  Energy  (K.E.),  Potential  Energy  (P.E.) 
and  the  Total  Energy  (T.E.)  for  the  experiments 
N3216  (dash-cross  lines),  F3216(5)  (dash  lines), 

fir 

<.  P3216(2-3)  (solid  lines)  and  P3216  (dotted  lines), 

as  function  of  time  . 
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Fig.  12.  The  difference  between  the  mean  zonal  flow  U at 
about  50.6*N,  on  day  9,  apd  30  m/sec  for  differ- 

r-  * 

ent  experiments. 
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interval  between  application  of  the  filter.  While  the  diagram  indicates 
similar  difference  in  the  zonally  averaged  u-field  for  intervals  of  time 
corresponding  to  points  A and  B,  the  details  of  the  simulation  are, 
in  general,  different.  A more  important  conclusion  from  Fig,  12  that, 
it  shows  that  differences  produced  by  the  periodic  filter  are  very 
sensitive  to  the  frequency  of  applying  that  filter. 

From  the  above  results  we  form  the  following  picture.  At 
the  outset  of  the  flow  evolution,  the  troughs  and  ridges  tilt  to  transport 
momentum  from  middle  latitudes  towards  the  north  and  south.  At  the 
same  time,  through  non-linearities , energy  is  fed  to  smaller  scale 
waves.  The  energy  flow  to  the  small  waves  reaches  a maximum  then 
starts  to  flow  back  to  the  larger  waves  as  the  trough-ridge  tilt  reverses 
direction  to  give  an  opposite  transport  of  momentum  indicative  of  the 
stable  nature  of  the  flow.  The  application  of  the  filter  interferes  with 
this  process.  When  the  filtering  is  applied  very  often,  very  little 
energy  leaks  out  of  the  system  and  so  the  flow  is  forced  to  (almost) 
conserve  the  energy  in  the  larger  scales  and  thus  to  oscillate.  If  the 
filtering  is  applied  rarely,  the  smaller  scales  can  become  saturated 
and  initiate  the  reverse  flow  of  energy  which  leads  to  the  oscillation. 
However,  when  the  filter  is  applied  at  an  intermediate  frequency,  a 
significant  amount  of  energy  is  removed  from  the  system  at  the  small 
scales,  and  the  large  scales  will  tend  to  feed  still  more  energy  to  these 
scales  much  as  though  the  integration  were  beeing  restarted  after  each 


filtering . 


Since  the  process  of  an  oscillating  tilt  must  be  related  to 
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Fig.  13.  Same  as  Fig.  8 except  for  the  experiment 
P32l6(2-3).  with  east-west  filtering  (solid 
line)  and  north-south  filtering  (dashed  line) 
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a redistribution  of  energy  among  north-south  waves  we  carried  out 
two  further  experiments  where  filtering  was  applied  either  in  the  north- 
south  direction  or  the  east-west  direction  exclusively.  Fig.l3(A  and 
B)  gives  the  results  of  these  experiments.  Note  that  the  application 
of  the  filter  in  the  east-west  direction  only  permits  the  appropriate 
oscillation  in  the  tilt  of  the  phase  of  wavenumber  four  whereas  the 
application  of  north-south  filter  causes  the  persistent  tilt  that  has  been 
seen  in  previous  experiments.  There  are  some  other  effects  going  on 
in  these  experiments,  but  the  primary  effect  is  clear. 

We  are  led  to  conclude  that  the  scales  of  motion  filtered 
out,  in  low  resolution  experiments  with  periodic  filter,  play  an  impor- 
tant active  role  in  the  total  evolution  of  the  flow.  Further,  the  process 
requires  a flow  of  energy  from  the  smaller  scales  to  the  large  scales 
a possibility  which  cannot  be  envisaged  by  a filter  which  serves  as  an 
energy  sink.  We  will  take  up  this  point  later  when  we  discuss  the  result 
of  an  integration  of  an  unstable  initial  condition,  namely,  Haurwitz 
wavenumber  six. 


I V . 3 . 4 Overall  control  of  instability 

We  have  looked  into  the  effect  of  filtering  and  the  inclusion 
of  the  friction  term  on  the  detailed  behaviour  of  the  flow.  Here  we 
present  an  overall  view  of  how  well  the  different  schemes  control  the 
development  of  instability  which  is  the  main  reason  for  including  frict- 
ion or  filtering.  Table  3 indicates  to  the  nearest  day  how  long  a part- 
icular integration  proceeded  before  numerical  instability  developed. 
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We  note,  from  Table  3,  that  the  use  of  the  3-gridiength  filter  extends 
the  calculation  as  far  as  our  longest  run  without  indication  of  any 

instability.  In  the  experiments  with  a friction  term,  three  values  of 

. s' 

kinematic  coefficient  of  viscosity  were  tried;  >)l=10  , v2=l0  and 

1 4 x * 

V3=10  m/sec.  In  the  low  resolution  experiments,  the  low  value  of 


TABLE  3. 


Experiment  notation 

No . of  days 

before 

"blow-up" 

Experiment  notation 

No . of  days 

before 

"blow-up" 

N 3216 

16 

N 6432 

8 

F3216(4) 

17 

F6432(4) 

9 

F3216(5) 

18 

F6432(5) 

12 

F3216(6) 

>20 

P3216(2- 3),  P6432 
p6432(2-3)  and 
P3216 

>20 

friction  managed  to  delay  the  numerical  instability  for  less  than  one 
day.  The  medium  value  of  friction  was  sufficient  to  add  another  day 
of  simulation.  The  high  value  of  friction  controlled  the  instability  to 
the  20  day  limit  of  integration,  but  produced  far  too  much  smoothing. 

* In  middle  latitudes  for  high  resolution  experiments  and  with  mean 
flow  of  30  m/sec  >)l,  ^2  and  >)3  give  values  for  the  local  Reynold's 

number  of  the  order  10*,  lo'  and  10  respectively. 


N. LATITUDE  (degrees)  N. LATITUDE  (degrees) 
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In  fact,  after  10  days  with  high  friction  coefficient  the  amplitude  of 
wavenumber  four  for  the  V-field  at  about  28,lPN  ( the  maximum 
value  ) decreases  to  about  47%  of  its  initial  value  and  to  63%  alter 
15  days,  while  without  friction  the  same  point  decreases  by  less  than 
2%  after  15  days  for  P3216  and  less  than  0.002%  for  P6432. 

Fig.  14(A,B,C  and  D)  shows  the  u-field  on  the  days  3 
and  6 for  the  experiments  N6432  and  F6432(5).  The  two  models 
are  very  similar  up  to  3 days  without  indication  of  instability.  On 
day  6 irregularities  are  apparent  in  low  latitudes  in  experiment 
N6432.  In  fact,  some  indication  of  irregularities  could  be  seen  as 
early  as  day  4.  These  irregularities  grow  rapidly  with  time  and 
cause  a "blow-up"  on  day  8.  The  experiment  N6432  was  repeated 
with  half  the  time  step  ( &t=2  minutes),  but  the  same  instability 
appeared  at  the  same  time  ending  the  calculation  on  the  same  day. 
F6432(5),  on  the  other  hand,  has  no  indication  of  instability  on  day  6. 
It  is  clear,  apart  from  the  areas  of  irregularities,  that  F6432(5) 
and  N6432  are  in  good  agreement  cn  day  6.  Fig.  5(D)  shows  the 
u-field  on  day  9 for  the  experiment  F6432(5).  The  appearance  of  a 
similar  irregularities  is  shown  in  the  u-field  of  that  day.  The  gro- 
wth of  these  irregularities  causes  the  "blow-up"  of  calculation  on 
day  12. 

The  instability  of  the  low  resolution  experiments  appears 
to  be  different.  The  first  indication  of  trouble  can  be  noticed  as  early 
as  day  10  at  which  time  significant  amplitude  is  contained  in  zonal 
wavenumber  two  at  all  latitudes.  This  behaviour  is  strange  since 
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TIME  (days) 


TIME  (days) 


Fig.  15.  The  logarithmic  distribution  of  the  amplitude  of 
wavenumber  two  as  function  of  time  at  the  latit- 
udes 5. 6°N  (•),  39. 4°  N (x)  and  84. -fN  (o).  The 
upper  and  lower  curves  are  for  the  experiments 
N3216  and  F3216(5)  respectively. 
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neither  the  interaction  of  the  primary  wave  (zonal  wavenumber  four) 
nor  aliasing  should  develop  zonal  wavenumber  two.  Fig.  15  shows 
the  distribution  of  the  amplitude  for  zonal  wavenumber  two  as  funct- 
ion of  time  for  the  experiments  N3216  and  F3216(5).  The  irregula- 
rities grow  with  time  and  finally  end  the  calculation  on  days  16  and 
18  for  N3216  and  F3216(5)  respectively.  Note  that  the  other  expe- 
riments, including  N6432  and  F6432(5)  , shows  practically  zero 
magnitude  for  wavenumber  two  at  all  latitudes  as  far  as  the  integra- 
tion goes . 

IV.  4 Experiments  with  wavenumber  6 

In  section  IV.  3 a stable  wave  was  studied  and  we  saw 
that  because  of  the  reversable  cascade  of  energy  and  enstrophy  betw- 
een large  and  small  scale  waves,  the  periodic  filter  with  low  resolu- 
tion produced  rather  large  differences  in  the  evolution  of  the  flow. 

In  this  section,  we  perform  similar  studies  with  initial  conditions 
drawn  from  a Haurwitz  wavenumber  6.  As  shown  by  many  authors 
this  wave  breaks  down  in  a few  days  forming  cut  off  lows  in  middle 
latitudes  (see  for  example  Merilees  (1974)  and  Kasahara  (1977)).  With 
this  initial  condition  it  is  unwise  to  run  the  experiments  with  3-grid- 
length  filter  for  low  resolution  (11  = 32  and  JJ=16),  the  initial  disturbance 
has  considerable  energy  in  wavenumber  12.  Merilees  (1974)  showed 
that  due  to  the  unstable  nature  of  the  flow  and  the  resulting  cascade 
of  enstrophy  to  the  shorter  wavelengths,  experiments  with  resolution 
up  to  15  waves  are  not  enough  to  predict  the  evaluation  of  this  wave 
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— N6432 
*«*  F 6432 (5) 
•••  P6432 (2-3) 
P6432 


ZONAL  WIND  (m/sec) 


MERIDIONAL  WIND  (m/sec) 


Fig.  16.  (A);  The  zonal  average  of  the  u-field  after  3 days  as  a function 

of  latitude.  (B ) : The  amplitude  of  wavenumber  6 for  the 
V-field  after  3 days  as  a function  of  latitude.  Thin  dash 
line  is  the  initial  profile. 
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after  3 days.  In  this  section  the  experiments  N6432,  F6432i3), 
P6432(2-3)  and  P6432  will  be  discussed.  In  addition  a control  run 
with  filtering  each  time  step  and  11=128  and  JJ  = 64  was  performed.  The 
control  run  will  be  considered  as  the  true  solution.  The  time  step 
used  in  these  experiments  is  four  minutes  except  for  the  control  run 
where  the  time  step  is  5 minutes. 

In  each  of  the  above  five  experiments  the  disturbance  soon 
developed  a north-west/south-east  tilt  which  transported  angular 
momentum  out  of  the  middle  latitudes.  The  five  models  are  in  essen- 
tial agreement  up  to  3 days  with  eastward  phase  progression  of  about 
23.7°  per  day  near  the  equator  and  21.6°  per  day  in  high  latitudes 
compared  to  24,6  per  day  predicted  by  the  non-divergent  barotropic 
theory.  Fig.  16(A  and  B)  shows,  on  day  3,  wavenumber  6 for  the 
V-field  and  the  mean  zonal  flow  for  the  experiments  N6432,  F6432(5), 
P6432(2-3)  and  P6432.  After  3 days  the  differences  between  the  above 
mentioned  models  become  more  noticeable.  All  the  models  indicate 
more  rapid  motion  of  the  waves  in  low  latitudes  after  4 days.  The 
main  difference  is  in  middle  latitudes  where  experiments  with  Fourier 
filter  show  a stationary  or  weak  westward  motion  of  the  waves  at  these 
latitudes  after  four  days.  The  other  experiments  indicate  weak  east- 
ward motion.  This  is  shown  on  Fig.l7(A  and  B)  where  we  plot  the 
distribution  of  the  phase  angle,  of  the  primary  wave  for  the  V-field, 
as  function  of  time  at  the  latitudes  30.9°N,  53. 4* N and  70.3*N.  As 
a result  of  the  difference  in  phase  speed  with  latitude  between  exper- 
iments with  Fourier  filter  and  the  other  experiments,  the  tilt  of  the 
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the  wave  (which  leads  to  the  export  of  westerly  momentum  in  middle 
latitudes)  is  steeper  in  P6432(2-3)  and  P6432  than  F6432(5).  This 
result  in  weaker  easterlies  for  F6432(5)  which  can  be  seen  on 
Fig. 18(A).  The  latitudinal  distribution  of  the  mean  zonal  flow  and  the 
main  wave  for  the  V-field,  on  day  6 ^Fig.l8(A  and  B)  and  Fig.l9(A 
and  B)3  , indicate  that  the  experiments  with  Fourier  filter,  on  that 
day,  are  closer  to  the  control  run  than  the  other  experiments.  N6432 
shows  clear  irregularities  on  the  mean  zonal  flow  and  the  main 
disturbance  of  the  V-field,  Fig.l8(A  and  B).  N6432  shows  an  increase 
of  kinetic  energy  of  about  28%  on  day  6 and  subsequently  "blows -up", 
Fig.  20.  F6432(5)  shows  irregularities  in  low  and  high  latitudes  on 
day  8,  Fig.  21(A),  and  the  calculation  is  ended  after  9 days  due  to  the 
accumulation  of  energy  in  high  wavenumbers.  At  this  time  the  kinetic 
energy  has  increased  by  about  53%,  as  shown  in  Fig.  20.  Fig.  21 
shows  the  h-fieids,  on  day  8,  for  the  experiments  F6432(5),  P6432(2-3), 
P6432  and  the  control  run.  Apart  from  the  area  of  irregularities 
shown  on  the  field  of  F6432(5),  the  comparison  between  the  fields 
shows  better  result  for  the  runs  with  the  Fourier  filter.  For  example, 
in  middle  latitudes  F6432(5)  shows  open  waves  while,  as  expected 
from  the  results  of  others  and  also  from  our  control  run  (Fig.  21(D)), 
P6432(2-3)  and  P6432  give  closed  systems  in  middle  latitudes. 

It  is  clear  from  this  section  that  the  periodic  filter  does 
control  the  aliasing  instability  and  at  the  same  time  gives  reasonable 
evolution  of  the  flow  with  time  with  little  additional  computing  time. 

On  the  other  hand  for  wavenumber  6 the  friction  term  with  coefficient 
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of  lO^m’/sec  is  not  enough  to  prevent  the  accumulation  of  energy  in 
high  wavenumbers  and  higher  values  of  the  friction  coefficient  will 
cause  excessive  smoothing. 

From  these  experiments  we  cannot  conclude  that  a non- 
linear viscous  term  would  not  do  a better  job,  but  certainly  the 
results  of  Miyakoda  et  al  (1971)  indicate  that  such  an  approach  is  not 
too  promising. 
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Fig.  20.  The  percentage  change  from  initial  time  of  the 
Kinetic  Energy  (K.E.)  for  the  experiments 

N6432  ( ),  F6432(5)  ( ),  P6432(2-3) 

( - ....  ),  P6432  (xwxx  ) and  the  control  run 
( • • • • ),  as  function  of  time . 
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Fig.  21.  The  free  surface  height,  afterS  days,  for  the  models 
F6432(5),  P6432(2-3),  P6432  aend  the  control  run;  A, 

B,  C and  D respectively.  Contour  interval  is  240  meter 


- 67  - 


CHAPTER  V 


SUMMARY  AND  CONCLUSIONS 


One  of  the  problems  of  long  range  forecast  models  is  their 
need  for  some  sort  of  smoothing  to  prevent  the  non-linear  (aliasing) 
instability  and  to  achieve  reasonably  smooth  fields  . Two  methods  are 
tried  for  a one  level  primitive  global  model.  The  use  of  a Fourier 
filter  which  eliminates  the  3-gridlength  waves  and  smaller  from  time 
to  time,  and  the  addition  of  a linear  friction  term,  in  the  momentum 
equation,  to  dissipate  the  kinetic  energy  especially  from  small  scale 
waves.  The  models  are  tested  with  two  types  of  simple  initial 
conditions,  a stable  wave  where  the  fields  tend  to  return  back  to  their 
initial  values  and  unstable  wave  where  the  fields  break  down  within 
few  days.  Up  to  three  or  four  days  the  models  show  no  significant 
differences  between  the  two  approac'ns.  In  fact,  even  runs  without 
any  smoothing  are  in  good  agreement  with  them. 

The  3-gridlength  filter  completely  controls  any  instability 
up  to  the  20  day  limit  of  our  runs.  It  also  requires  little  additional 
computer  time.  The  time  evolution  of  the  flow  is  quite  reasonable 
except  if  the  filter  interferes  with  an  essential  exchange  of  energy 
between  the  shorter  and  longer  waves.  In  this  case  the  filter  can 
cause  serious  deviations  from  the  true  solution.  In  this  case  the 
solution  to  the  error  is  to  increase  the  resolution  of  the  model. 

On  the  other  hand,  a coefficient  of  viscosity  equal  10  mVsec. 
fails  to  prevent  the  non-linear  instability  while  a value  of  10* m*/sec 


smoothed  the  fields  too  much.  Further,  the  use  of  the  friction  term 
requires  considerably  more  computer  time  as  indicated  in  section  II.  8. 
Perhaps  a nonlinear  viscosity  formulation  would  be  more 
useful,  but  based  on  the  work  of  Miyakoda  et  al  (1971)  this  does 
not  appear  promising. 

Somewhat  surprisingly,  the  application  of  the  3-gridlength 
filter  each  time  step  gave  the  best  overall  results  for  both  the  stable 
and  unstable  wave.  Such  a procedure  is  however  wasteful  of  resolution 
when  the  pseudospectral  algorithm  is  used  for  the  estimation  of 
derivatives.  Further,  such  a procedui  > would  effectively  convert 
the  pseudospectral  model  into  a spectral  model  with  about  2/3  of  the 
number  of  degrees  of  freedom.  In  that  case  it  may  well  be  more 
efficient  to  use  a spectral  model  from  the  start. 

While  the  results  generally  indicate  that  numerical 
filtering  is  better  than  a friction  term,  it  is  not  obvious  that  the 
particular  filter  used  is  best.  In  fact,  some  less  drastic  form  of 
numerical  smoothing  may  be  able  to  provide  equally  good  long  term 
numerical  stability  without  strongly  interfering  with  energy  exchanges 
in  a stable  type  flow. 

Finally,  it  should  be  noted  that  these  results  have  been 
obtained  with  very  simple  initial  conditions.  As  such,  they  can  only 
be  indicative  of  sort  of  conclusion  that  would  be  obtained  with  initial 
conditions  from  real  data. 
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Similarly  we  can  prove  that 
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Appendix  (II) : CFL  criterion 


In  order  to  obtain  guidelines  for  the  choise  of  time  step 
we  treat  the  numerical  stability  problem  in  Cartesian  rather  than 
spherical  coordinates.  This  approach  was  taken  by  Jacques(1976) 
with  useful  results.  Here  we  include  a friction  term.  Consider  in 
Cartesian  coordinates  the  following  equations 


and 


4.  — ^ , 


where  ^ =gH  is  a constant. 

Using  Fourier  series,  let  us  express  all  the  dependent  variables  in 
the  form 

A - 111  A0  6.  > 

A is  function  of  time,  k is  the  east-west  wavenumber  and  i is  the 
north-south  wavenumber.  Using  a leapfrog  time  scheme  and  reme- 
mbering that  the  space  derivatives  are  exact,  the  finite  differences 
analogues  to  the  above  equation,  for  a particular  mode  k and  1 , 
can  be  written  as 


t+6t  t-fit 

U.  - 

^*1  <•* 


t -1  at  [ ik  4?  + \)  ] , 

_ 1 it  [ i 1 <£.' 

-ii  <|>  C k.  -t- 1 n//  3 
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where  we  have  used  the  standard  technique  of  calculating  the  friction 
term  at  time  t-  &t. 

Assuming  solutions  of  the  form 


•w 

O- 


= ^ 


lult 

e . 


where  u,  V?  and  4>*  are  constants,  we  find  that 

v*(a-a'J  -v  3i)  • 

c^CA-a')+2«  At# -H* 

iU  at 

where  A = 

The  condition  for  a non-trival  solution  is  that  the  determinant  of  the 
algebraic  system  vanish,  i.e. 

C.A-A'  ) + 2^  A C k'+  £.  ) O 

© (A->‘)+tv>AtA’^v+r; 

«*.*  fct  ^ k x i At«£  i 

If  we  denote  A = (A“A)‘*‘^^  <ltA  (^T^)> 

then  expansion  of  the  determinant  yields  the  condition 
A | ( A -A  ) A + 4 (At)1-  £ (kl  +l')^  =0. 

In  order  for  the  calculation  to  be  stable  we  must  have  all  possible 


j/.Atk 

X\AXi 


U-A) 


=o 


* Case  1. 

The  first  case  we  consider  is  when  the  first  factor  in  the  above 
condition  is  zero,  i.e. 

( A ) + 2 $ At  ( k*+  )b'  =0, 

or 

= 1 - 2$  At  ( k*-+£). 

Then  to  insure  stability  we  require  that 

-1  1 - 2$  At  ( k*  + JLl  ) < 1. 

The  r.h.s.  is  always  satisfied.  The  l.h.s.  gives  a condition  for 

[ : 

stability 

At  ^ 1/  $ ( k*+*0. 

t C 

Without  friction  these  modes  are  neutral,  with  friction  we  expect 
the  modes  to  be  exponentially  damped.  Note  however  that  if 

l/^k'+T)  At  ^ 1/2$  (kVio, 


then  £ < 0 and  the  solution  is  not  physically  realistic  since 
there  will  be  a sign  change  each  second  time  step.  Thus  a reason- 
able condition  to  impose  is 

At  ^ 1/2$  (k‘+J.‘). 

Case  2 . 

When  the  second  factor  is  zero  we  have 

( X - X f + 2 $ At  ( k*  + r ) a'  ( A - A1  ) + 4(  At)\§  ( k*  +£V)  = 0, 
or  (AX-lf+2$  At(k*+il)  (jf-i)  + 4 (Atff  ( k*  +!' )A*  *0. 

i 


1 
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Putting 

a = \>  &t  ( kfc+  1*  ) and 

the  above  equation  can  be  written  as 


b = 2(  Mt  § ( k*  + f ), 


( A1  -1  )*■  + 2a  ( A1  -1  ) + 2b  {jf  -1  ) + 2b  = 0 , 


or 

X*  + 2 ( « 

where 

X 

ii 

So 

X = - 

or 

r— * 

II 

Suppose 

that 

a r b 


— | * 

) - 2b] 


( a + b ) - 2b  C 0 , 

then  the  stability  condition  j \ ^ \ implies  that 


[(  l * a - b f -{(a+bf  - 2b^] 


1 , 


or 


k 


( 1 - 2a  j1  ^ 1 . 


By  definition  a ^ 0 and  since  we  have  choosen  a ^ 1/2  , from 
the  study  of  the  first  factor,  we  find  that  this  condition  is  always 
satisfied.  But  if 

a + b ) - 2b  ^ 0 , 


t 


l 
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will  contribute  to  decreasing  the  middle  expression.  It  is  therefore 
necessary  and  sufficient  that 

1 - (a  + b ) - [(a  + b f - 2b^t’>  - 1 , 


or 


( a + b ) + 


D 


X 

) 


a + b ) - 2b 


-A 

U* 


It  follows  that 

r x ^ 

|_  ( a + b ) - 2b  J ^ 2 - ( a + b ) . 

Now  if  ( a + b ) > 2 , then  the  inequality  can  never  be  satisfied, 
therefore 


2 - ( a + b ) must  be  0 . 

Then  it  is  possible  to  square  the  previous  inequality,  which  gives 


2 ^ 2a  + b , 

i £ 0 ( at ) ( kl  + X1  ) + ( a t ) # ( k*  + iv ) . 

Without  friction  this  condition  reads 

At  1/  ( kl  + £ , (36) 

V, 

where  C = ( g H ) *•  is  the  speed  of  gravity  waves. 

In  our  experiments  H = 8 km  and  the  largest  value  used  for  $ , 
the  kinematic  viscosity,  is  10*  m/sec.  Also  the  smallest  value 
used  for  At  , in  our  experiments,  is  one  minute.  This  means  that 
in  our  experiments  At  <£>  is  greater  or  of  the  same  order  as  ^ . 
So  with  fully  explicit  scheme  we  can  say,  in  our  experiments,  that 
the  friction  do  not  have  significant  effect  on  the  CFL  criterion  which 
can  be  considered  (36). 


r 


Let  us  now  integrate  the  same  system  of  equations  using 
the  one  dimensional  semi-implicit  scheme  discussed  in  II. 4.  Using 
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a leapfrog  time  scheme,  then  for  a particular  mode  k and  1 we  get 

-2  &tQk  ^ d ( k*  £*+  jTu^  )”], 

V^-  ■\£tV*  = -2  &tQ$  4>*  +0  ( k*-  + $.'  Q , 

&*  4*=-zi  At§[k^%^  + ukl. 

Assuming  solution  in  the  form  (35),  the  equations  become,  with 

. itdfct 

-A  = e , 

ufA-A*  )+i  atk^fA+A  )tZ^  at  A*  (kx+i  )u*=0, 

■/(A  -A*  ) + 2i  At  5.  4?  + 2>)  at  a*  ( k*  + £l ) v*  = o , 

<£(A  -A*  ) + i at§  k ( a + a'  ) I?  + 2i  at  £ "**  = 0 . 


u* 


V*  and  4^  have  non-trivial  solution  if 


( A -a’  ) + 2S)  at  A ( kl  + l*  ) 


.-I 


( A -A1  )+2  ^ atA*  ( k*-  + i1 ) 


i at  k (A  +A  ) 


2i  at£ 


r\ 


i # at  k ( A +A  ) 


2i  at§-i 


(A  -A  ) 


= 0, 


With  A = ( A -A*  ) + 2 >)  at  j[  ( k1  + 5.  ),  we  get 

a[(A  -a'  ) a+H(  at  + (at*  k1^  ( A +a'  )]  = o. 

Which  gives 

A = 0, 

with  the  condition  of  stability,  as  discussed  before,  namely 

at  1 / 2 0 ( k‘  + £ ) . 


The  other  solution  gives 


Mi  ■ 
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( A - A*  * + 2 ) ( At  )X  ( A - X ) ( k'  + ^ ) + 4 ( A t ) §-  !*■  + 

+ ( ) k1-^  ( A +A*  )■  = 0 , 


°r  ( £ -l  ) + 2 $ ( At  ) ( k%  +?  ) ( A*  - 1 ) + ( At  ) k1-  § ( A + l S’  + 

+ 4 ( At  ?<§•**’ A*  = o . 

Putting 

a = ^ at  ( k*  + £ ) , b = ( l\  t ) k*«£  and  c = 2 ( & t ) § 5.  , 

the  previous  equation  reads 

X (l+b)  + 2X(a  + 2b  + c)+2(2b+c)  = 0, 


where 


or 


x = A - 1 • 

X — — ( a.  + 2b  + c ) _+  [“  + 2b  + c ) - 2(l+b)(2b+cT[,t|/(Ub) 
£ = 1 + j T \ Wy  | - (a+2b+c)+  (a+2b+c?  -2(l+b)(2b+c)J4^ 


Following  the  same  procedure  as  before,  we  get  as  a condition  for 
stability 

1 >.  0 ( At  ) ( kl  +!'’)  + ( At  *§>8.  . 

The  condition  without  friction  is 

At  1 / He  . 


- 82  - 


IBfVP'V  ■ 


Appendix  (III) : The  relation  between  the  momentum  transport 

and  the  tilt  of  the  V-field 


In  the  northern  hemisphere  a north  (south)  transport  of 
westerly  momentum  is  associated  with  north-west/south-east  (north- 
east/south-west) tilt  of  the  V-field.  This  relationship  between  the 
tilt  of  the  V-component  of  wind  and  the  direction  of  momentum  tran- 
sport is  strictly  valid  only  in  the  case  of  gecstrophic  winds.  Here 
we  perform  a sample  calculation  of  the  actual  momentum  transport 
to  indicate  that  this  relationship  is  valid  for  the  numerical  experim- 
ents described. 

The  momentum  transport  across  a given  latitude  circle 
is  proportional  to  ^ 


\ 

iT 


^ -V  dA 


> 


where  u and  V are  the  deviations  of  the  velocity  components  from 
the  zonal  average.  In  the  sample  calculation,  by  far  the  dominant 
component  is  wavenumber  four.  Therefore  we  can  write 


u = A cos  (4  A - S>* ) , 

V7  = B cos  (4  A - f ^ ) , 

where  A and  B are  the  amplitudes  for  \J  and  V*.  \Sf^  and  ^ 
are  the  phase  angles. 

It  follows  that  the  momentum  transport  can  be  written  as 


NORTHWARD  TRANSPORT  I SOUTHWARD  TRANSPORT 


0 90 

PHASE  DIFFERENCE  (degrees) 


Fig.  22.  The  latitudinal  distribution  of  the  difference  between 
the  phase  angles,  for  wavenumber  4,  of  the  u and  V 
fields  for  the  experiment  F3216(5).  After  3 days 
(solid  line),  6 days  (dash  line)  and  9 days  (dot  line). 


N. LATITUDE 


— - uV  dA  = — — \ AB  cos(4>  - v?  ) cos(4X  - ^ ) dA 

2.TT  ) J 

» ° 

AR 

= — Y~  cos(^-^  ) . 

The  previous  equation  shows  that  the  direction  of  the  momentum 
transport  across  a latitude  circle  is  determined  by  the  phase  differ- 
ence between  u'  and  V with  the  following  rules: 

Northward  transport  of  westerly  momentum  for 

/ 2 < ^ -W/  2 , 

with  maximum  at  *4  = O • 

Southward  transport  of  westerly  momentum  for 

TT  / 2 < V?*  - ^ < 3T  / 2 , 

with  maximum  at  =TT  . 

The  momentum  transport  is  zero  for,  ="^*/2  or  3*^/2. 

Fig. 22  shows  the  distribution  of  with  latitude 

on  days  3,  6 and  9 for  the  experiment  F3216(5).  Comparing 
Fig. 22  with  Fig.3(A)(  page  33  ),  it  is  clear  that  we  get  the  proper  sign 
of  the  momentum  transport  by  inspecting  the  tilt  of  the  v-field. 


